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In experiments the two-dimensional systems are realized mainly on solid substrates which intro¬ 
duce quenched disorder due to some inherent defects. The defects of substrates influence the melting 
scenario of the systems and have to be taken into account in the interpretation of the experimental 
results. We present the results of the molecular dynamics simulations of the two dimensional sys¬ 
tem with the core-softened potential in which a small fraction of the particles is pinned, inducing 
quenched disorder.The potentials of this type are widely used for the qualitative description of the 
systems with the water-like anomalies. In our previous publications it was shown that the system 
demonstrates an anomalous melting scenario: at low densities the system melts through two continu¬ 
ous transition in accordance with the Kosterlitz-Thouless-Halperin-Nelson-Young (KTHNY) theory 
with the intermediate hexatic phase, while at high densities the conventional first order melting 
transition takes place. We find that the well-known disorder-induced widening of the hexatic phase 
occurs at low densities, while at high density part of the phase diagram random pinning transforms 
the first-order melting into two transitions: the continuous KTHNY-like solid-hexatic transition and 
first-order hexatic-isotropic liquid transition. 

PACS numbers: 61.20.Gy, 61.20.Ne, 64.60.Kw 


Despite almost forty years of investigations, the con¬ 
troversy about the microscopic nature of melting in two 
dimensions {2D) still lasts. In a crystal in contrast to an 
isotropic liquid two symmetries are broken: translational 
and rotational. These two symmetries are not indepen¬ 
dent, since a rotation of one part of an ideal crystal with 
respect to another part disrupts not only the orienta¬ 
tional order but also the translational order. However, it 
is possible to imagine the state of matter with orienta¬ 
tional order, but without the translational one. As it was 
shown by Mermin [ij in two dimensions the long-range 
translational order can not exist because of the thermal 
fluctuations and transforms to the quasi-long-range one. 
On the other hand, the real long range orientational order 
does exist in this case. 

These ideas were used in the widely accepted KTHNY 
theory [H where it was proposed that, in contrast to 
the 3D case where melting is the first-order transition, 
the 2D melting can occur through two continuous transi¬ 
tions. In the course of the first transition the bound dis¬ 
location pairs dissociate at some temperature Tm trans¬ 
forming the quasi-long range translational order into the 
short-range one, and long-range orientational order into 
the quasi-long range order. At this transition the decay of 
the translational correlation function will change from al¬ 
gebraic to exponential, and the orientational correlation 
function will obey the algebraic decrease. The new in¬ 
termediate phase with the quasi-long range orientational 
order is called the hexatic phase. After consequent dis¬ 
sociation of the disclination pairs at some temperature 
Ti the hexatic phase transforms into the isotropic liquid 
with the exponential decay of the orientational correla¬ 
tions. 


The KTHNY theory seems universal and applicable to 
all system, however, despite numerous experimental and 
simulation studies, it could not be consistently verified or 
refuted. The KTNHY scenario was unambiguously ex¬ 
perimentally confirmed for superparamagnetic colloidal 
oarticles interacting via long-range dipolar interaction 
On the other hand, in 2D the first-order melt¬ 
ing without hexatic phase is also possible [Il|) lUl • The 
possible mechanisms of the first-order melting are the 
formation of the grain boundaries 11| and the disso¬ 
ciation of the disclination quadrupoles 12|. It is now 
widely believed that the melting process strongly de¬ 
pends on the pair interaction of the system. Numerous 
experimental and simulation studies demonstrate that 
the systems with very short range or hard core poten¬ 
tials melt through weak first-order transition, while the 
melting scenarios for the soft repulsive particles favor the 
KTHNY theory. However, controversies still exist even 
for the same systems, like, for example, for hard spheres 


Recently, a number of papers have appeared where 
another melting scenario was proposed [2J-l30|. It was 
argued that the basic hard disk model demonstrates 
the two-stage melting transition with a continuous solid- 
hexatic transition but a first-order hexatic-liquid transi¬ 
tion 24 - 2^ . In Ref. 271 it was shown that in 2D Yukawa 


system there is a two-phase coexistence region between 
the stable hexatic phase and isotropic liquid. Kapfer and 
Krauth have analyzed the effect of the potential soft¬ 
ness on the melting scenarios. They have proposed that 
the soft sphere system with the potential of the form 
U{r) = {ojr)'^ melts in accordance with the KTHNY 
theory for n < 6, while for n > 6 the two-stage melting 
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transition takes place with the continuous solid-hexatic 
transition and the first-order hexatic-liquid one. 

In experiment, planar confinement is typically realized 
by adsorption on solid substrates which introduce frozen- 
in (quenched) disorder due to some roughness or intrin¬ 
sic defects. Quenched disorder can influence the crys¬ 
tallization/melting scenario of two-dimensional systems. 
As it was shown in Refs. 31, SH (see also 33, 3^), the 
KTHNY melting scenario persists in the presence of weak 
disorder. The temperature of the hexatic-isotropic liquid 
transition Ti is almost unaffected by disorder, while the 
melting temperature Trp , drastically decreases with in¬ 
creasing disorder 3ll434l|. and the stability range of the 
hexatic phase widens. Recently these predictions have 
been tested by experiment and simulations of the behav¬ 
ior of the superparamagnetic colloidal particles 

In this paper, we present a computer simulation study 
of 2D phase diagram of the previously introduced core 


softened potential system [35l - l4fll | in the presence of 
the quenched disorder. The core-softened potentials 
are widely used for the qualitative description of sys¬ 
tems demonstrating the water-like anomalous behavior, 
including density, structural and diffusion anomalies, 
liquid-liquid phase transitions, glass transitions, melting 


maxima [351-14 

The system we study in t he p resent simulations is de¬ 
scribed by the potential 35-4^ 47-4^: 


!/(?') =e(^) -f (1 - tanh(fci{r - CTi})). (1) 

where n = 14 and kiu = 10.0. a is the hard-core diame¬ 
ter. We simulate the systems with the soft-core diameter 
tJi/cr = 1.35. As it was shown in our previous publica¬ 
tions , the phase diagram of the system consists of 

three different crystal phases (see also similar picture in 
0), one of them has square symmetry and the other two 
are the low density and high density triangular lattices 
(see, for example. Fig. 6 (a) in [i^). Melting of the low 
density triangular phase is a continuous two-stage transi¬ 
tion, with a narrow intermediate hexatic phase, in accor¬ 
dance with the KTHNY scenario. At high density part of 
the phase diagram, the square and triangular phases melt 
through one first order transition. At higher tempera¬ 
tures and high density there is one first order transition 
between trianglar solid and isotropic liquid. 

In the remainder of this paper we use the dimensionless 
quantities, which in 2D have the form: f = r/cr, P = 
Pa^je^ V = VjNa^ = l/p,T = ksT/e^di = cti/ct. In 
the rest of the article the tildes will be omitted. 

We use the molecular dynamics simulations of the 
system in NVT and NVE ensembles (LAMMPS pack¬ 
age [HI) with the number of particles equal to 20000. 
Quenched disorder is introduced by pinning a randomly 
chosen subset of particles to the random positions and 
let them to be immobile for the entire simulation run. 
We make the simulations of 10 independent replicas of 


the system with different distributions of random pinned 
pattern and after that average the thermodynamic func¬ 
tions over replicas. We calculate the pressure P versus 
density p along the isotherms, and the correlation func¬ 
tions G(,{r) and Grir) of the bond orientational ■i/'e and 
translational ipr order parameters (OPs), which charac¬ 
terize the overall orientational and translational order. 

We define the translational order parameter 'ipT 
(TOP), the orientational order parameter ipQ (OOP), 
the bond-orientational Ge{r) (OCF) and translational 
Grir}^ (TCF) correlation functions in the conventional 
way 


Q (TCF) correlation functic 


TOP can be written in the form 



( 2 ) 


where is the position vector of particle i and G is 
the reciprocal-lattice vector of the first shell of the crys¬ 
tal lattice. The translational correlation function can be 
computed in accordance with the definition Si,®: 


Grir) 


< exp(iG(rj 


gir) 



(3) 


where r = — r^l and g{r) =< 6{ri)S{rj) > is the 

pair distribution function. The second angular brackets 
< ... >rp correspond to the averaging over the random 
pinning. In the solid phase the long range behavior of 
Grir) has the form [1,|J Grir) oc r~^'^ with \ < px < \- 
To identify the existence the orientational order and 
the hexatic phase, we define the local order parameter, 
which measures the 6-fold orientational ordering, in the 
following way: 


^'6(l’i) 


1 

n{i) 



1=1 


(4) 


where 9ij is the angle of the vector between particles i 
and j with respect to a reference axis and the sum over 
j is counting n(i) nearest-neighbors of j, found from the 
Voronoi construction. The global OOP is obtained as an 
average over all particles and random pinning: 



(5) 


The orientational correlation function Ge{r) (OCF) is 
given by the equation analogous to Eq. jS) 




( 6 ) 


where 'k6(r) is the local bond-orientational order param¬ 
eter (H]). In the hexatic phase there is a quasi-long range 
order with the algebraic decay of the orientational corre¬ 
lation function GQ{r) oc with 0 < pe < 3 M- 
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FIG. 1: (Color online) (From top to bottom) Translational 
correlation functions (a) and enveloping curves (b) for: 1). no 
random pinning; 2). concentration of pinning centers equal 
to 0.1%; 3). 0.2%; 4). 0.5% (p = 0.48, T = 0.15). 


Before to proceed with the study of the phase diagram 
of system o in the presence of the random pinning, let 
us consider the influence of disorder on the orientational 
and translational correlation functions. In accordance 
with the results of Nelson and coworkers 31, ^ 


did not find any qualitative change in the behavior of 
Ge{r) (Eq. ([ 6 ])) in the presence of pinning. In contrast 
with Ge{r), we found the qualitatively different behav¬ 
ior of the translational correlation function Grir) (Eq. 
(IH)). In Eig. [T]we present Grir) for p = 0.48, T = 0.15 
without pinning and for 3 different values of the con¬ 
centration of pinning centers: 0.1%, 0.2%, 0.5% (Concen¬ 
tration 0 . 1 % corresponds to 20 particles frozen-in in the 
random positions in the system). Density and temper¬ 
ature p = 0.48, T = 0.15 correspond to the point on 
the phase diagram deeply inside in the low density solid 
phase. Without pinning we have the conventional power 
law decay of TCF. In the case of pinning, the slope of 
the enveloping line increases at some crossover value tq. 
The region r < vq corresponds to the local order unaf¬ 
fected by random impurities, while asymptotic behavior 
of TCE for r > ro is determined by the random pinning 


31j . From Fig.[I](b) one can see that, in accordance with 


the intuitive picture, rp decreases with increase of impu¬ 
rities concentration along with the increase of the slope of 
the enveloping line. Further we will consider the equal¬ 
ity Vt = 1/3 for the long range asymptote of TCF as 
the solid-hexatic transition criterium. The hexatic-liquid 
transition occurs at 775 = 1/4. 


Let us first consider the effect of pinning on the low 
density part of the phase diagram where melting occurs 
in accordance with the KTHNY theory. Fig. [2] illustrates 
the behavior of OCF and TCF for different temperatures 
at p = 0.48 and the concentration of the pinning centers 
equal to 0.1%. Using equations pr = 1/3 and pe = 1/4, 
we calculated the low-density phase diagram in the pres¬ 
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FIG. 2: (Golor online) OCF and TGF for different temper¬ 
atures at p = 0.48 and concentration of the pinning centers 
0 . 1 % 


ence of random pinning. The phase diagram is shown in 
Fig. |3] (a). One can see that random pinning leaves al¬ 
most unaffected the line of the liquid-hexatic transition, 
while considerably changes the location of hexatic-solid 
transition. As a result, in the presence of random pin¬ 
ning the hexatic phase drastically widens in accordance 
with the theoretical predictions 31, 321 (see also iE3)- 


The effect of random pinning on the first-order melt¬ 
ing transition is shown in Fig. [3] (b). The temperatures 
higher than 0.3 are considered. In this temperature range 
the transition in the pure system with the potential o 
is of first order |47 h 49I | . The lines determined from the 
equations px = 1/3 and pe = 1/4 are denoted by the 
open symbols. One can see, that the line where GQ{r) 
decays algebraically with the exponent equal to 1 /4 is lo¬ 
cated inside the two-phase region of the first-order melt¬ 
ing transition of the system without pinning. At the same 
time, in the presence of random pinning the line of solid- 
hexatic transition obtained from condition px = 1/3 is 
shifted to higher densities. From this picture one can 
conclude that random pinning transforms the single first- 
order transition into two transition with rather wide hex¬ 
atic phase. The solid-hexatic transition is continuous, 
while the nature of the transition from the hexatic phase 
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FIG. 3: (Color online) Low density (a) and high density (b) 
phase diagram of the system with the potential © withont 
pinning (filled symbols) and with the concentration of the 
pinning centers equal to 0.1% (open symbols). 


to the isotropic liquid is unclear. In order to elucidate 
this issue, we plot the isotherms of the system in Fig. |4] 
with and without random pinning. One can see, that 
the van der Waals loops which are characteristic of the 
first order transition are almost unchanged by impuri¬ 
ties. This means that the hexatic phase transforms into 
isotropic liquid through first-order transition. 

In conclusion, in this paper we present the computer 
simulation study of melting transition in 2D core soft¬ 
ened system in the presence of random pinning. It is 
shown that at low densities, where the system without 
pinning melts through two continuous transitions in ac¬ 
cordance with the KTHNY scenario, random disorder 
does not change the character of the transition but dras¬ 
tically widens the range of the hexatic phase. At the 
same time, at high densities where the conventional first 
order transition takes place without random pinning, dis¬ 
order drastically change the melting scenario. The single 
first-order transition transforms into two transitions, one 
of them (solid-hexatic) is the continuous KTHNY-like 
transition, while the hexatic to isotropic liquid transition 
probably occurs as the first order transition in accordance 


with the recently proposed scenarios [24 - 


The pos- 


FIG. 4: (Color online) Isotherms of the system with the po¬ 
tential © without pinning (open symbols) and with the con¬ 
centration of the pinning centers equal to 0.1% (filled sym¬ 
bols) for T = 0.3 and T = 0.4. 


sible mechanism for this transition is the spontaneous 
proliferation of grain boundaries ESIiS- 

It should be noted, that the nature of the first-order 
liquid-hexatic transition is still puzzling. As it was shown 
in BE[13, using the Landau-type expansion one can 


see that the solid-liquid transition can be of first order in 
the case of small enough thermal fluctuations and con¬ 
tinuous KT-like, when the singular fluctuations of the 
phase of the order parameter takes place. The choice be¬ 
tween these two possibilities depends on the form of the 
potential and the thermodynamic parameters. In con¬ 
trast, in the case of the liquid-hexatic transition one can 
show that both Landau expansion and KT mechanism 
lead to the continuous transition, however, simulations 
predict the first order transition. The strict theory like 
the KTHNY one is necessary for explaining this contro¬ 
versy. The results of this study can be useful for the 
qualitative understanding the behavior of water confined 
in the hydrophobic slit nanopores 52, 5^. 
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